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ABSTRACT 


Anode sheaths impact the operation of many practical plasma devices. This 
complex region is explored in detail for collisional, isothermal (identical specie 
temperatures), low-temperature plasmas, where sheath dimensions are in the micron 
range. The selected approach involves postulation of a specific electric field 
distribution with two shape factors. Previous research regarding planar anodes is 
verified and expanded upon using greater parameter ranges. ‘z', a dimensionless 
quantity specifying plasma composition and condition, groups diverse plasmas into 
‘families’ exhibiting similar sheath characteristics. 'n', a nondimensional ratio of 
electrical energy to thermal energy in the sheath, allows temperature effects to be 
studied. The investigation focuses on three disparate plasma families that span a z 
range of 1.1729 to 2.1493, at n values defined by plasma temperatures of 6000°K, 
3000°K, and 300°K. Results indicate that at lower temperatures, charge production 
in the outer sheath is generic to the electric field distribution, and that the sheaths 
themselves are nearly unaffected by substantial changes in temperature (1.e., n). 
Conversely, sheath density and extent are shown to vary significantly for differing z 
values. Newly-derived equations governing cylindrical anodes generate sheaths that 
are virtually identical to corresponding planar cases. It is shown that only those 
anodes whose radii are comparable to the plasma's ‘characteristic radius’ (y) must be 
treated with the cylindrical formulation; non-vacuous plasmas would require micron 
width anodes to be thus affected. Finally, an analytical approach yields solutions that 
confirm the numerical results, and offers an algebraic approximation for high-n 


plasmas. 
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Il. INTRODUCTION 


Various types of plasma devices have been in operation for over 30 years. 
Marshall's coaxial plasma gun [Ref. 1] was successfully accelerating volumes of 
hydrogen plasma as far back as the late 1950's. In the intervening years, the number 
of applications for plasma devices has steadily increased to include 
magnetohydrodynamic (MHD) power generation, laser pumping, strategic defense, 
and electromagnetic propulsion for interplanetary spacecraft. However, due to their 
somewhat complex nature, comparatively little is understood about the sometimes- 
destructive sheath regions surrounding the electrodes in every plasma device. This 
work attempts to further previous research efforts concerning description of the 
anode sheath. 

As stated by Biblarz [Ref. 2], completely satisfactory anode sheath solutions do 
not exist for several plasma conditions: one such case involves steady, collisional, 
low-temperature, isothermal discharges. He then goes on to derive an involved, 
nonlinear differential equation that describes the entire plasma region affected by a 
planar anode, from the surface to the undisturbed plasma. A presumed (but 
‘shapeable') function describing the electric field in that region, selected after much 
deliberation [Ref. 2 and Ref. 3], ultimately allows charge production rates and 
electron/ion populations to be plotted as a function of distance from the anode. 
Nondimensional parameters make the solution profiles applicable to several families 
of plasmas. 

In this work, numerical techniques are employed to verify the one previously 
solved case, and to explore planar anode solutions to several other plasma conditions. 
This can be readily accomplished due to the generality of the formulation and the ease 


with which the profiles can be produced. Sheaths and ambipolar regions are properly 


generated for all cases. The effect of nondimensional parameter variations on the 
sheath is extensively investigated, and some general conclusions are hypothesized. In 
addition, an analytical technique is presented that supports the numerical results and 
allows for an accurate algebraic solution to low-temperature conditions. Finally, 
similar derivations produce equations that treat the cylindrical anode sheath problem; 
the equations and their profiles are then analyzed, and a comparison is made to planar 


anode findings. 


Il. THE CARTESIAN PROBLEM 


This section discusses and validates one previous anode sheath research effort 
concerning steady, low-temperature collisional plasmas. In addition, several new 
parameter cases are examined and analyzed. The geometry of the Cartesian problem 


is illustrated in Figure 1. 


V 


(distance from the surface) 


Figure 1 Cartesian Anode Geometry 
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A. DERIVATION REVIEW 

Previous work by Biblarz [Ref. 2] treats the planar anode sheath problem in a 
one-dimensional Cartesian fashion, and is reviewed here as it forms the starting point 
for this work. Electric field intensity (E), electron and ion population densities ( n, 
and n, ), and all other relevant quantities are considered to vary only with linear 
distance (y) from the flat anode surface. The applicable relations are Gauss’ equation 


and the two species continuity equations, as shown on the following page. 


(Ee fo-nd (1a) 


dy 
dn. 
i -eunE-eb{ 2] (1b) 
dy 
dn 
Vea Cllennied ep{ 2. (1c) 
dy 


Note that e is the electron charge constant, €, is the permittivity constant, the 
subscripts e and 1 indicate electron and ion particles, the j's are the species 
contributions to the total current density (i.e., J=j,+j,), the u's are the particle 
mobilities, and the D's are the diffusion coefficients. Although this set does not fully 
constrain the problem (as an equation describing plasma reactivity, n,, is absent), any 
selection of a specific form for E(y) implicitly fixes that variable.. It is of interest to 
explore the problem with ‘guesstimates’ of E(y) in order to facilitate solutions. 

Reference 2 then introduces two new parameters, designated K° and K , which 
are defined by Equations (2a), (2b), and (2c). These K-terms are somewhat artificial 
parameters, although they are related to the total current. More meaningful is the 
derivative of K’, which is directly proportional to the net production rate of 


charges, n, (i.e., GK) =n_/D, ). 


ea = 238 (2a) 
eb. eD 
. ea (2b) 
ED, cP) 
2J 
and K’=#-K +— (2c) 
eD 


Manipulation of the previous six equations and the Einstein relation (which 
relates the mobilities to the diffusion coefficients) yields the isothermal differential 
equation below. Note that k is the Boltzmann constant, and that all primes denote 


derivatives with respect to the variable y. 


kh KR es 08 p {Ms se (3) 
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Advanced knowledge of E(y) considerably facilitates the solution of Equation (3). 
Biblarz' justification for the selected form of E(y) is detailed elsewhere [Ref. 2 and 
Ref. 3]. This function meets several critical conditions inherent to the stated 
problem, including the required distribution of E and the behavior of both K-terms at 


the boundaries. The function is reproduced as Equation (4) below. 


| A 
E(y) =E, Pl Tay (4) 


The parameters ‘a’ and 'A' are shape factors, where ‘a’ is of the order of the 
sheath length, and ‘A' relates to the physical parameters Ad, and E, (as discussed on 
pages 7 and 8). 

Substitution of the chosen E(y) into Equation (3), followed by an order-of- 
magnitude analysis, the nondimensionalization of the parameters in accordance with 
Figure 2, and further simplification based on field and current properties at the 
electrode and in the undisturbed plasma, all combine to yield the governing 


differential equation and boundary condition of Figure 3. 
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Figure 2 Nondimensionalized Parameters 





Figure 3 Simplified Nondimensional Differential Equation 
and Boundary Condition Governing the 1-D Cartesian Problem 


The introduction of the dimensionless parameters n and z further generalizes the 
problem, in that one solution to the above equation can be applicable to a large 
number of specific dimensional cases. Their additional significance is discussed in a 
subsequent section. 


As will be seen, both numerical and analytical solutions to the equations of Figure 


3 are possible, yielding profiles of K” and (K") that are a function of distance from 


the planar anode. 


Further manipulation of Equations (1) and (2) derives relations that produce 
electron and ion population profiles [Ref. 2]. A more precise form of these equations 
is given in Figure 4. Note that specific solutions for K* are required to generate the 


n-profiles. Several distinct cases in the next section yield population curves that 


clearly illustrate the propriety of the sheath and ambipolar regions in the plasma near 
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the anode. 
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Figure 4 Ion and Electron Population Equations 





Reference 2 also presents two more useful relations from analysis involving 
physical observations (Equations (5) and (6)). Recent analysis has proven that the 
infinite series in Equation (6b) converges for all values of z. These are combined and 
manipulated to yield the two important equations of Figure 5, which allow 
determination of the parameters z and n. The procedure 1s as follows: 


¢ the type of gas defines the anode potential drop Ad, (which is essentially 
equivalent to the ionization potential of the gas) 


¢ the particular case or application specifies E,, and n.. 


¢ zis then found by iteration of the implicit equation at the top of Figure 5 


¢ the choice of T, locks in n using the bottom equation in Figure 5 


a {7a exp(2z’ ) (5) 
en_ 
Ad, =E.VA f(z) (6a) 
oc qe 
where f(z) =z 2fam-i)mi (6b) 
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Figure 5 Equations That Yield z and n 





Specific values of z and n allow computation of the electric field parameters A 


and a (Equations (5) and (6)), which in turn yield the tailored form of the E 


distribution (Equation 4), and ultimately permit the generation of K’, (K*) , and n.,, 
profiles in the plasma regions close to the anode (equations of Figure 3 and Figure 4). 
B. NUMERICALLY SOLVED CASES 
1. Procedure 
Numerical solution of the governing differential equation (Figure 3) can be 


achieved with a fourth-order Runge-Kutta FORTRAN program. First the dependent 


variable is redefined and the equations are rewritten in accordance with Figure 6. 


<[w] = 2n- {E hy eau 








dy E, (¥+1) 
w(Q) = 
where w(y) = Fo 


Figure 6 Modified 1-D Cartesian Differential Equation 
and Boundary Condition for Numerical Solution 


Application of the Runge-Kutta scheme computes the value of w at each y. 
The data for K* is then recovered with the division of each w datapoint by the 
corresponding value for E(y)/E,. (K*) data are extracted with the following steps 


for each value of y: 
* compute W using the defining differential equation (Figure 6) 


¢ perform the operation below (derived from the Quotient Rule) 


drs) adwl &) ox E) Ee 
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Finally, the n,; profiles are computed as a function of y using the available 


& 


information and the equations of Figure 4. A single computer routine can be made to 
perform all of the required operations (see Appendix A); the resulting data are then 
plotted. 

Note again that this technique is substantially simpler than alternate methods 
that do not presume a specific form of E(y). Ensuing discussions address the 
advantages and disadvantages of this approach. 

2. Case I: the Nitrogen Problem 

The first specific case involves verification of an earlier example [Ref. 2]. A 
planar anode in contact with nitrogen plasma was analyzed under the following 
conditions (nitrogen is a common discharge gas with well-known properties): 

* nitrogen's anode potential drop (A@,): 15.51 V (singly ionized) 

¢ the electric field strength in the undisturbed plasma (E,.): 12,000. V/m 

¢ the particle density (n,): 10°° m™ 
¢ the temperature: 6000°K 


The expression used to approximate f(z) (Equation 6a) was truncated after 13 
terms, producing an error in the calculation of z which is on the order of 10°%. 
Using the given data, the following parameters result: 

7 = iSOZ6 

i) =e 

© a= 1.95421x10° m and A =1.17793 x10" m° 
¢ E, = 262,265. V/m 

The computed value for E, (obtained using Equation 4) relates specifically to 
all plasmas that are defined by the above values of z and n; this includes nitrogen 
plasma at the stated conditions. The normalized electric field E(y) for this case is 
shown in Figure 7. Note that the field decreases monotonically and abruptly to a 


constant value in the undisturbed plasma, as 1s required of the model. 


/ 


The corresponding K’, (K*) , and n,, profiles are presented in Figure 8 and 
Figure 9. Note that those Figures also represent all other plasma cases whose gas 
composition, density, temperature, and electric field intensity are defined by 
Z=1.75626 and n=99.12285. Figure 9 clearly illustrates the distinct sheath and 
ambipolar regions present for this general case.” The K-term profiles of Figure 8 


likewise exhibit expected characteristics, with K’ rising roughly monotonically from 


+1 to +2, while (K*) decreases to near-zero at the outer edge of the sheath. The 


leftmost downturn of (K*) , while possibly correct, is quantitatively suspect since the 


model of ionization by electron impact breaks down at the fringe of the collisionless 


region (i.e., for y/a < ~10°). 


* The sheath is that small region, extending away from the anode, where the 
electron and ion populations are not equivalent; in the ambipolar region, the ratio 
of oppositely-charged particles is 1:1, but the populations of both are less than those 
that exist in the ‘undisturbed plasma’ 
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Figure 7 Normalized Electric Field (z=1.75626, n=99.1229) 


11 


K-term Profiles (Case I) 
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Figure 8 K-term profiles (z=1.75626, n=99.1229) 
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n-Profiles (Case I) 
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Figure 9 n,, profiles (z=1.75626, n=99.1229) 
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As the parameter ‘a’ specific to Case I is known, it can be seen that the 
nitrogen plasma sheath for the particularly defined discharge extends approximately 
0.1 mm out from the planar anode. Similarly, the region of undisturbed plasma for 
such a device begins approximately 6 mm from the anode surface. Both of these 
results are in accordance with expectations, the generalities of which were obtained 
after much deliberation [Ref. 4]. 

These numenically-obtained Case I results confirm the work of Reference 2. 
Postulation of a suitable ionization/recombination mechanism is left to that discourse. 


The remaining objectives of this work are as follows: 


* investigation of the effects which are produced as the defining parameters are 
varied over a wide yet practical range 


* presentation of an analytical solution to the general Cartesian anode sheath 
problem 


¢ derivation and examination of the related Cylindrical anode sheath problem 
3. Cases II and III: Varying Temperature in the Nitrogen Problem 
Comprehensive investigation of the nitrogen plasma problem requires 
consideration of temperature's effect on the anode sheath. Thus, while all other Case 
I conditions are held constant, two somewhat more practical values for the isothermal 
temperature are considered: 
¢ Case Il fixes T, at 3000°K 
- Case III sets T, at300°K (this is representative of discharges in laser pumping) 
Temperature-only variations impact n exclusively; z and the normalized 
electric field E(y)/E, (including the inherent factors A and a) are unaffected. Since 
n is by definition inversely proportional to T, (Figure 2), the decreased temperatures 
of Cases If and HI signify larger values of n. As a result, the first-order derivative 
term in the governing differential equation (Figure 3) becomes less influential at the 


lower temperatures. Indeed, for the conditions of Case III, the contribution of the 


14 


derivative term is negligible (note, however, that this circumstance does not nullify 
the validity of the formulation; the equation does not become degenerate). 

Application of the lower equation in Figure 5 yields n values of 198.246 and 
1982.46 for Case II and Case III, respectively. The K-term profiles for these general 
cases are presented in Figures 10 and 11. Both cases exhibit the same general K-term 
tendencies that have been recognized previously, with minor variations that illustrate 
the decreasing importance of the net-ionization term. 

The effect of temperature variation (or, equivalently, n variation) on a 


plasma in the anode region is explored using Figures 12 and 13, which directly 
compare the K’ and (K’) profiles of the three cases. One significant observation, 


contributed by Biblarz subsequent to Reference 2, hypothesizes that the asymptotic 
behavior of the K-terms with decreasing temperature reflects the approach to 
ionization that is generic to the electric field distribution; the overall results are thus 
independent of the details of ionization and recombination. As an important practical 
example, measurements of charge density would not reflect relevant details of the 
ionization/recombination mechanisms. 

Conversely, the n., profiles for this general z=1.75626 case are nearly 
unaffected by the changes in temperature. Figure 14 illustrates charged particle 
curves that exhibit only a barely perceptible shift of position over the entire range of 
temperatures tested. This unexpected result has substantial implications; the most 
significant of these is the acceptability of the isothermal-particle assumption (1.e., the 
assertion that T. = T, = T,). Though such a premise may be unrealistic, its validity 
appears to be inconsequential to the population profiles. Electric-field effects 
dominate any temperature effects in the cases examined. In fact, the physical 
significance of n (stated in Reference 2 as a ratio of electrical to thermal energy) 


predicts such dominance with its large magnitudes. 
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K-term Profiles (Case II) 
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Figure 10 K-term profiles (z=1.75626, n=198.246) 
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K-term Profiles (Case III) 
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Figure 11  K-term profiles (z=1.75626, n=1982.46) 
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Kplus Profiles for Varying Temperatures 
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Figure 12 K’ profiles (z=1.75626) 
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(Kplus) Profiles for Varying Temperature 





Figure 13. (K°) profiles (z=1.75626) 
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n-Profiles for Varying Temperature 
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Figure 14 n., profiles (z=1.75626) 
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C. EXPLORING THE ENVELOPE FOR z AND n 
The foregoing results include a preliminary investigation into the effect of n- 
variation on the planar anode sheath. This upcoming section more fully explores how 
changes in the parameters z and n impact the sheath region. 
1. Practical Extremes of 7 
The parameter z is a dimensionless representation of three basic variables 
(Ad , E. , and n,,) that specify the condition of a plasma; several different gases, at 
differing densities and field strengths, can be related by their equivalent z values. 
Various diverse plasmas belonging to the same 'z-family’, at the same temperature, 
would theoretically exhibit identical sheath characteristics. Using the technique 
discussed previously, anode sheaths for plasmas defined by two limiting values of z 
are examined. 
Singly-ionized elements possess anode potential drops (A@,) that range from 
24.46 Volts (helium) down to 3.87 Volts (cesium), with the low values being 
preferred. Common densities (n,) for these type collisional plasmas vary 
approximately from 10'* to 10” particles per cubic meter. Finally, reasonable field 
strengths in the undisturbed plasma (E,,) are somewhat arbitrarily set from 5000 to 
30,000 Volts per meter. Application of the upper equation in Figure 5 yields the 


following ‘practical’ extremes for z: 


Pigg s 2 2.1495 (8) 

This relatively small numerical range represents an extremely diverse 
grouping of plasmas; z is a decidedly sensitive parameter. Small z values represent 
those plasmas with small potential drops and low densities that are subjected to very 
high electric fields (e.g., low density cesium at 30,000. V/m). Plasmas with large z 


values are characterized by gases with large potential drops, at conditions of high 


2) 


density and low field strength (e.g., high pressure helium at 5000. V/m). By way of 
comparison, the nitrogen plasma considered earlier in this work (median potential 
drop and density, small field strength) has a slightly-larger-than-median z value of 
1.7563. 

Values for n are required to completely constrain the problem, which entails 
the choosing of one or more temperatures. To facilitate direct comparisons, the same 
arbitrary quantities for 1, that have been considered for the nitrogen case are 
designated as the standard values (i.e., 6000°K, 3000°K, and 300°K). 

Consequently, numerical sheath solutions are next generated and examined 


for these forementioned limiting conditions: 
¢ ‘Small z' plasmas (z= 1.1729) at the three standard temperatures (n= 16.4135, 
32.8269, and 328.269 for T, at 6000°K, 3000°K, and 300°K) 


¢ 'Large z' plasmas (z= 2.1493) at the standard temperatures (n= 257.563, 
515.125, and 5151.25 as T, decreases to 300°K) 


2. Sheath Solutions for ‘Small z' Plasmas 
K-term profiles for the 'small z' plasmas (at the three standard temperatures) 
are depicted in Figures 15, 16, and 17. Qualitatively, each of these curves displays 
the same previously noted characteristics inherent in the nitrogen plasma curves: 


monotonically increasing and well-behaved K*, and a rate-of-charge-production term 


(K*) that decreases and vanishes at the outer edge of the sheath. As before, the 


leftmost downturn of (K"*) iS quantitatively suspect at the fringe of the collisionless 
region. 

These reassuring results extend to the composite curves of Figures 18 and 19, 
which depict K-term profile variations as a function of temperature (1.e., n) changes. 
In particular, the charge production rate curves of Figure 19 again illustrate the 


striking implication that charge production in the outer sheath is independent of the 
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K-term Profiles, Low z, T=6000 K 
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Figure 15 K-term profiles (z=1.1729, n=16.4135) 
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K-term Profiles, Low z, T=3000 K 
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Figure 16 K-term profiles (z=1.1729, n=32.8269) 
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K-term Profiles, Low z, T=300°K 
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Figure 17 K-term profiles (z=1.1729, n=328.269) 
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Kplus Profiles, Low z, Various Temp. 
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Figure 18 K’° profiles (z=1.1729) 
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(Kplus) Profiles, Low z, Various Temp. 
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Figure 19 (K*) profiles (z=1.1729) 
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temperature. Moreover. the inner sheath once again sees an increase in the net 
production of charges as the temperature decreases. 

Electron and ion population profiles for the ‘small z' plasmas are given in 
Figure 20. As with the nitrogen case (and other plasmas of that z-family), the sheath 
and ambipolar regions are clearly visible and quite appropriate, while the population 
curves themselves again appear to be nearly unaffected by large variations in 
temperature. 

Such qualitative observations for the ‘small z' plasmas indicate that previous 
impressions concerning ‘median z’ plasmas (e.g., nitrogen) are not limited to that one 
category, but are possibly characteristic of all collisional low-temperature plasmas. 
The sheath profiles generated for the ‘large z' cases of the next section further 
confirm this hypothesis. 

Quantitative deviations in the sheath due to the change of z are addressed 1n a 
later section. 

3. Sheath Solutions for 'Large z' Plasmas 

The high density, small field ‘large z' plasmas continue the trends established 
by the other collisional plasmas, albeit with some minor differences. One such 
difference is illustrated by the K-term profiles of Figures 21, 22, and 23. The curves 
display the same general traits of the other plasmas, with the exception that the K* 
curve no longer climbs monotonically. This small anomaly becomes more 
pronounced at higher temperatures (Figure 21), yet does not appear to significantly 
impact the charge production rate curves, which still behave as expected. 


Similarly, there are no surprises in the composite curves of Figures 24 and 


25. In fact, the (K’) curves of Figure 25 demonstrate even more strikingly the 


generic, electric-field-dependent nature of ionization and recombination in the sheath. 


28 


n-Profiles, Low z, Varying Temp. 
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Figure 20 n,, profiles (z=1.1729) 
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K-term Profiles, High z, T=6000°K 
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Figure 21 K-term profiles (z=2.1493, n=257.563) 
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K-term Profiles, High z, T=3000° K 
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Figure 22 K-term profiles (z=2.1493, n=515.125) 
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K-term Profiles, High z, T=300°K 
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Figure 23 K-term profiles (z=2.1493, n=5151.25) 
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Kplus Profiles, High z, Various Temp. 
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Figure 24K’ profiles (z=2.1493) 
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(Kplus) Profiles, High z, Various Temp. 





Figure 25 (K*) profiles (z=2.1493) 
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The n-profiles of this class of plasmas are depicted in Figure 26. While there 
are some yet-to-be-addressed differences from previous n-profiles, the sheath and 
ambipolar regions are again evident, as are the profiles’ theorized independence from 
temperature changes. 

4. The Effect of z on the Planar Anode Sheath 

Results to this point have revealed some important and perhaps non-intuitive 
concepts concerning the properties of planar anode sheaths, concepts which appear 
applicable to all low-temperature, collisional plasmas. However, sheath variations 
attributable solely to changes in the parameter z have yet to be discussed. 

Comparison of the n-profiles for each z-family (Figures 14, 20, and 26) 
illustrates a subtle point: although all n-profiles have been shown to be nearly 
independent of temperature, they do become slightly more sensitive to temperature 
variations as the value of z decreases. This is possibly attributable, in part, to the 
lower densities of ‘small z’ plasmas; the temperature-dependent kinetic energy 
changes of the particles may produce a more observable effect on the total 
ionization/recombination process in a plasma not already saturated with density- 
driven collisions. Perhaps more contributory is the previously-discussed dominance 


of electrical energy over thermal effects. The ‘small z' plasmas exhibit smaller 


magnitudes of E, than other plasma families, decreasing such dominance. The 
involvement of Ad, in this phenomenon remains unclear. In any event, such n 
profile temperature shifts remain nearly negligible for all plasmas considered to this 
point, which further validates the innocuous nature of the constant temperature 
assumption. 

Figure 27 depicts the K” curves of three diverse z-families at the same 


temperature (6000°K). As noted previously, all are well-behaved and do not differ to 
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n-Profiles, High z, Varying Temp. 
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Figure 26 n., profiles (z=2.1493) 
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Kplus Profiles, 6000°K, Various z 
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Figure 27 K’* at6000°K for Different z-Families 
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any significant degree. Contrastingly, the charge production rate curves of these 
three plasma families (Figure 28) are somewhat disparate and very insightful. ‘Large 
z plasmas show significant charge production throughout the entire sheath, in 
contrast to the localized and decreased production rates depicted for the ‘small z’ 
plasmas. This may indicate the fact that electric-field effects are more localized for 
smaller z. 

The n-profiles for the three z-families are depicted in Figure 29. Unlike 
similar profiles for varying temperature (i.e., n), these curves are decidedly affected 
by changes in z. Most prominently, the magnitude of the charged particles changes 
dramatically as z varies, especially in the sheath. This is a graphic indication of z- 
induced changes in both E, and current flow for diverse plasmas. Also noteworthy 1s 
the fact that the sheath itself extends further from the anode surface as the value of z 
decreases. In particular, the sheath of the ‘small z' plasmas stretches over ten times 
the distance occupied by the ‘large z’ sheath. 

5. Some Final Thoughts on the Influence of n 

The effect of n variations on the anode sheath has already been indirectly 
addressed via the extensive consideration of temperature's influence on the problem. 
The bottom equation of Figure 5 inversely relates these two parameters for fixed z 
and A@d,. Physical interpretation of n as a ratio of electrical and thermal energy has 
also been reviewed. 

The purely numerical importance of n is apparent with a glance at the 
goveming differential equation (Figure 3); its reciprocal controls the influence of the 
charge production rate term. Large values of n reduce this term to near zero, the 
consequences of which can be seen in every low-temperature case (i.e., 300°K). To 
fully and completely explore this parameter's effect on the sheath, solutions are 


presented for the original nitrogen plasma at artificially small values of n (in 
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Figure 28 (K*) at 6000°K for Different z-Families 
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n-Profiles as a Function of z, T=6000 K 
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particular, n is set equal to 25., 10., and 1.0). Note, however, that the validity of 
these curves is questionable at best, since such small n values equate to some 
extraordinarily high temperatures for the nitrogen (up to 595,000°K!), the existence 
of which contradicts the original assumption of low-temperature plasmas. 
Nevertheless, the following curves may perhaps reveal some valid trends. 


Despite the unusual and probably erroneous K* plots for the small-n 


conditions (Figure 30), the (K*) curves of Figure 31 hearteningly retain some of the 
familiar traits concerning charge production rate in the sheath. The corresponding n 
profiles are presented in Figure 32; remarkably, only the anomalous 'n=1' profiles 


exhibit tendencies somewhat contrary to those already observed. 


D. ANALYTICAL SOLUTIONS 

In order to confirm the validity of the numerically-obtained solutions, an 
analytical solution to the problem's governing relation (Figure 6) is presented. 

1. Procedure for 'The Outer Expansion Method' 

Reference 5 offers an approach which is ideally suited to the form of the 
differential equation in Figure 6. The term (1/n) is justly defined as a ‘small’ 
parameter, and the equation is already in dimensionless form. A series solution for w 
is presumed, approximate yet valid for ‘non-small’ values of y. Equation (9) depicts 


the power series expansion for w: 


k ie 3 
a l | l l 
w= Ywi-| =w,t—iw,ti—|wotl—iw,t: (9) 


Substitution of Equation (9) into the governing equation allows terms with 


like powers of (1/n) to be equated, which in turn allows each w, to be solved for 


analytically. Because the solution for w, is purely algebraic in this formulation, the 
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Figure 30 K’ profiles (z=1.75626, Small n‘) 
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(Kplusy’ Profiles, z for Nitrogen, Low eta’s 
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Figure 31 (K‘) profiles (z=1.75626, Small n°‘) 
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n-Profiles, z for Nitrogen, Low eta’s 
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first-order derivatives in each of the subsequent w, expressions can be evaluated 


exactly (i.e., symbolically). The total solution w is thus an infinite summation of 
exact expressions. Prevalent magnitudes of n allow highly precise approximations of 
w with the series truncated to only three terms. 


Utilization of this technique produces the following approximate analytical 


‘outside’ solution (w,) to the governing differential equation of Figure 6: 


sw.cwe{t| {4}. 
Ma a ae ae n W> 
where: vf 2.0 acca 
LEG VE.) (G+) 
~ [CEG WE) A Es (EGE M9 +I | Ltt) 
-| 6.0 |(22)/-| 2.0 [sy ; 12.0 r 
eiey | E, (E(y)/E,) h E, (E(y JE, )( ¥ +1) 


oe EO | (E(G)/E,). | 
(E(¥)/E,) (¥ #1) (ECH)/E My +1) 


























Figure 33. Three-Term Approximate Analytical Solution (w, ) 


The structure of the equations themselves offers insight. Note that each w, is 
a function only of z and y; they are completely independent of n (i.e., temperature). 


The contribution of n to the total solution manifests itself solely in the series- 


expanded expression for w,. Note also that the equation for w, is identical to the 


original governing differential equation for the case when n is allowed to increase to 


infinity. As a result. the partial solution w, is itself a valid approximation for the 
true solution (w) at large values of n. This fact is also corroborated by the series- 
expanded expression for w,,; all terms beside w, vanish for large n. 

As before, plots of K* can be recovered from the data for w (or, in this case, 
w,). Using parameter values from the first nitrogen plasma case (z=1.7563 and 
n=99.1229), K* profiles from both the ‘Outer Solution’ and the Fourth-Order Runge- 
Kutta scheme are compared (Figure 34). Although the ‘Outer Solution’ does indeed 
appear to diverge for extremely small values of y, the two curves are nearly 
identical. Such correlation offers comforting evidence for the validity of the 
numerical procedures and results presented in this work. 

Figure 35 illustrates the soundness of using just w, to recover data for large- 
n cases. The solid curve is from the Runge-Kutta solution for the 300°K nitrogen 
case (z=1.7563 and n=1982.5); the dotted profile depicts data generated using only 


the partial solution w,. 
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Kplus (Case I): Numerical vs. Analytical 
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Figure 34 Comparison of Numerical and Analytical Solutions 
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Kplus (High eta): Numerical vs. Analytical 
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Figure 35 One Term (w,) Algebraic Approximation for High n°‘ 
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Ii. THE CYLINDRICAL PROBLEM 


The following sections present detailed derivations and analysis for the 
cylindrical anode sheath problem. As with the planar cases, the focus 1s limited to 
steady, low-temperature collisional plasmas. Exploration of the cylindrical sheath is 
desirable, as several functioning plasma devices contain or employ a cylindrical 
architecture. In addition, these results can be compared and contrasted to the 
previously-presented Cartesian results. The geometry of the one-dimensional 


cylindrical problem is illustrated in Figure 36. 
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Figure 36 Cylindrical Anode Geometry 


A. DERIVATIONS 

Attainment of a governing differential equation for the cylindncal case begins 
with the re-derivation of the applicable relations (Gauss' equation and the two species 
continuity equations) in cylindrical form. Reference 6 presents the following 


differential form of Gauss’ law: 
Vapoee (10) 
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Vis the del operator, E is the vector form of the electric field, p is the charge 


density (which can be represented here in accordance with Equation 11), and e, is the 
permittivity of free space. The divergence (V-E) of the electric field in cylindrical 


coordinates is given in Equation (12) below: 


p= lnngeae) (11) 
gel) G 1\OE, dE, 
vB =( 7) 508) +(o ae is 


For a one-dimensional electric field that varies only as a function of radius (r), 


combination of Equations (10), (11), and (12) yields the following one-dimensional 


cylindrical form of Gauss' equation: 


dE (1 e 
—— m= —_—_ aia l 
= «(=e = mn) | (13) 


Comparison of Equations (13) and (1a) accentuates the appearance of a new term 
which 1s the result of the cylindrical derivation. 
The species continuity equations are derived from conservation equations [Ref. 7] 


that can be manipulated into the form given below: 


-j.=-epun.E-eD.Vn, (14a) 
j =eun.E -eD Vn, (14b) 


For a one-dimensional electric field, the previous relations are easily rewritten as 


shown in Equations (15a) and (15b) on the following page; these are the 1-D 
cylindrical species continuity equations. 
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d 
peng (15a) 


dn, 


j, =eun,E- eD, 
dr (15b) 


Note that these species equations are nearly identical to those for the planar case 
(Equations 1b and I|c); only the independent variable has been renamed. 


The procedure from this point is the same as that of Reference 2: 


¢ Combine these forms of Gauss' equation and the species continuity equations 
with the Einstein relation 


¢ Incorporate the K-terms as defined by Equation (2) 


* Assume isothermal conditions 


¢ Derive a single nonlinear differential equation in K’, E, and their derivatives 
The desired intermediate differential equation (the cylindrical counterpart to 
Equation 3) is presented below. Note that all primes here denote derivatives with 


respect to the variable r. 


kTfK\ . 2 fe Yi E 
SE} ao ee ee 
CHEE} E] « 
e° E E- rE rE rE fe 


In comparison to Equation (3), several new terms exist that are due solely to the 
cylindrical nature of the derivation. However, the validity of the above expression is 
demonstrated by letting r increase towards infinity, which causes the equation to 
approach the planar case. In that limit, Equation (16) reduces exactly to the 
corresponding Cartesian expression of Equation (3). 

The next step calls for an ‘educated guess’ of the form for E(r), and substitution 


of that function into Equation (16). For reasons stated previously, the basic form of 


a 


the electric-field function given in Equation (4) will be retained, although in a form 


modified for use with a cylindrical anode. The plasma electric field begins at the 


anode surface, which for a cylindrical anode is at r=r, (as opposed to r=0, the anode 


center). Thus, the presumed electric field function is 


B 
E(r) = E af S| (17) 


where 4 is the anode radius, B and b are constants that give specific shape to the 
field, and the range of ris from r, to ~. 
Following substitution of E(r) into Equation 16, the resulting lengthy expression 


is made dimensionless through use of the relations and definitions given in Figure 37. 
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Figure 37 Cylindrical Nondimensionalized Parameters 





Note that q is the cylindrical counterpart to z, and that r has a more conventional 
range of zero to infinity. The new parameter y appears because the simplifications 
that were previously employed (based on field and current properties at the electrode 
surface and in the undisturbed plasma) in the Cartesian derivation [Ref. 2] cannot 
eliminate all of the same constants when applied to the cylindrical case. y has units of 


length, and can be thought of as a ‘characteristic radius’ of the plasma conditions and 


anode width that define q and n.. 
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Further simplification of the differential equation is achieved with the derived 


J enE. 
eg KT, ) bat 


Finally, order-of-magnitude analysis reveals that some of the surviving terms 


expression [Ref. 2]: 


cannot be neglected unless the anode radius (r,) is sufficiently large. As a result, the 


differential equation can take one of two possible forms; the appropriateness of either 
is a function mainly of the anode's radial magnitude. Thus, the governing differential 
equation and boundary condition for ‘normal-sized' anodes (e.g., r,=10 mm) are 
given in Figure 38, while the corresponding equations for a ‘wire-thin' anode (e.g., 


r, ~0.1 mm) are given in Figure 39. 





Figure 38 Simplified Nondimensional Differential Equation 
and Boundary Condition, Governing the 1-D 'Normal Anode 
Cylindrical Problem 


t 





Figure 39 Simplified Nondimensional Differential Equation 
and Boundary Condition, Governing the 1-D 'Wire Anode' 
Cylindrical Problem 
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The equations of Figure 39 are actually appropriate for all magnitudes of 1, but 
some of the terms become negligible for the larger anodes. In fact, the degree of 
influence exerted by the extra terms of that differential equation (Figure 39) is 
directly tied to the ratio y/r, ; unless the cylindrical anode's radius is smaller than, or 
within an order-of-magnitude of, the ‘characteristic radius’ (y), then those extra 
terms are insignificant. For this reason, the simplified differential equation of Figure 
38 is also provided. 

Note also that this ‘normal anode’ differential equation is identical in form to the 
corresponding planar equation (Figure 3), which indicates that in spite of a tortuous 
derivation, most cylindrical and planar anodes disturb the plasma in nearly the same 
manner. Finally, as r in the expression of Figure 39 increases toward infinity 
(approaching the planar case), the extra terms vanish and the planar differential 
equation is recovered. 

The parameters q and n. have comparable meaning to their planar counterparts, 
and their specific values are computed using similar techniques. It can be shown 
through extensive manipulations that the cylindrical equivalents to Equations (5) and 


(6) are as follows: 


PAH : ; 

b -| 7h exp( 2q’) (19) 
Ao, =E. vB f(q) (20a) 
where f({q)= D> aa (20b) 
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Note that the infinite series converges for all values of q, and that the y/r, term 
in the denominator of Equation (19) is only relevant for the 'wire-thin’ (and smaller) 
anodes. The above equations are combined to generate the two important equations 
of Figure 40, which are used to compute the parameters q and n. in accordance with 
the following familiar procedure: 

¢ the gas composition defines Ad, 

¢ the particular case or application specifies E,. and n, 

¢ the device in use fixes the anode radius r, 

¢ q is then computed using the implicit equation at the top of Figure 40 


- the plasma temperature T, then yields n. via the bottom equation of Figure 40 
¢ for the extremely thin anodes, y is computed from its definition and the value of 


E(t = 0); both expressions are found in Figure 37. 


en,AQ, = 3 2»\_| Ad, f>,.2 
s.28. i q f(q)exp(2q°) Ae Jenn( 2a" 
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n.= eb Jew q f(q) 





Figure 40 Equations That Yield q and nm. 


As with the planar cases, specific values of q and n. allow computation of the 


electric field parameters B and b (Equations 19 and 20), which in turn yield the fitted 


form of the E distribution (Equation 17), and ultimately permit the creation of K’, 


aK ) , and n., profiles in the plasma regions close to the anode. All that remains to 


be derived are the cylindrical forms of the equations that compute those n-profiles, 
which are produced through manipulation of Equations (2), (13), and (15). The 


resulting elaborate expressions are presented in Figures 41 and 42. 
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Figure 41 Ton Population Equation 
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Figure 42 Electron Population Equation 


Gratifyingly, these n-profile equations produce their planar counterparts as 1, 
approaches infinity. Note also that the expressions for n, and n, are valid for all 
cylindrical anodes. Unfortunately, only one term in Figure 41 could be neglected for 
larger, 'normal-sized’' anodes, which stifles the notion of a separate ‘simplified’ 


equation for such cases. 


B. NUMERICALLY SOLVED CASES 
1. Procedure 

The previous derivations allow detailed numerical analysis of cylindrical 
anode sheaths in certain plasma conditions, using a procedure similar to that 
performed on the planar problem. As before, the governing equation is rewritten in 
terms of w, as defined in Equation (21). The resulting differential equations and 
boundary conditions for both the 'normal-sized’ anode (~5 mm radius and larger; this 
form also applies when (y/r,) << 1) and the 'wire-thin' anode (~0.5 mm radius and 


smaller, or when (y/r,)~1) cases are presented in Figures 43 and 44. Note that the 
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form of the equations for the 'normal-sized' anode 1s identical to that of the 


corresponding planar equations (Figure 6). 
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w(t) = E/E, (21) 
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Figure 43 Modified 1-D Cylindrical Equations 
for Numerical Solution (‘Normal-Sized' Anode) 





Figure 44. Modified 1-D Cylindrical Equations 
for Numerical Solution (‘Wire-Thin' Anode) 


As before, application of a FORTRAN Runge-Kutta algorithm (see Appendix 


A), combined with the previously-outlined data manipulations, yield the numerical 


profiles for K*, (K*), and n,, for any desired plasma and anode width. 
2. The Nitrogen Problem with a Cylindrical Anode 
The nitrogen plasma conditions designated previously as 'Case I' for the 
planar anode are now analyzed for both the 'normal-sized' and 'wire-thin’ cylindrical 


anodes. The Case I conditions are repeated on the next page: 


a7 


¢ nitrogen’s anode potential drop (A@,): 15.51 V (singly 1onized) 

¢ the electric field strength in the undisturbed plasma (E,): 12,000. V/m 
* the particle density (n..): 10°° m~ 

* the temperature (1,): 6000°K 

¢ the anode radii under consideration (r,): 10.0 mm and 0.1 mm 


These specifications and the newly-derived equations produce the following 


parameter values for their respective anodes: 


TABLE 1 CYLINDRICAL PARAMETERS FOR CASE I CONDITIONS 


ees Un lhe 
q— 1.75656 1.78774 
1.7 SMS 102.818 
v= N/A 3.9597=10° m 


BS 1.953x10% m 1.813x10°% m | 
Ee 262,543. V/m 293,216. V/m 


A companison reveals that the parameter values for the 'normal-sized' (10. 
mm) anode differ negligibly (only 0.1% and less) from their corresponding planar 
values. This is not an unexpected result, in light of the cylindrical equations’ 
similarity to their planar counterparts for the bigger anodes. More surprising are the 
values computed for the 'wire-thin' (0.1 mm) anode; while their differences from the 
planar values are not insignificant, only E, differs by more than 8%. The extreme 
sensitivity of the parameter gq, however, implies that percentage comparisons may be 


misleading. 
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Direct comparison of the planar and cylindrical profiles best illustrates the 
near-negligible differences in their respective sheath and ambipolar regions. Figures 
45 and 46 present the Case I K-term profiles produced by the three anode types 
(plane, ‘normal’ (10. mm) cylinder, and ‘wire’ (0.1 mm) cylinder), while Figure 47 
depicts the n-profiles for each of those anodes. As can be seen, profiles for the 
‘normal’ cylindrical anode are visually indistinguishable from the planar anode 
profiles. Amazingly, even the ‘wire’ anode profiles are nearly identical to those of 
the planar anode. This is a welcome result; despite many differences in the derivation 


and appearance of the governing equations for each type of anode, it seems that for 


this ‘small y' family of plasmas those cylindrical anodes whose radius is greater than 
that of a human hair can be treated to behave as planar. For such plasmas, only those 


cylindrical anodes whose radii are on the order of the sheath thickness itself would 
exhibit significant non-planar characteristics. 

To further validate this hypothesis, the three anodes are also compared under 
the lower-temperature ‘Case III’ circumstances, which represent the same conditions 
as Case I except that T, = 300°K. Consequently, the only parameter listed in Table | 


that changes is n.; its Case III values are 1983.14 for the 'normal-sized' cylinder, and 


/ 


2056.36 for the ‘wire’ anode. Figures 48, 49, and 50 show the K’, (K"), and n,, 


profiles (respectively) for Case III conditions. Once again, the planar and 10-mm- 
radius anodes disturb the plasma in identical fashion, while the 0.1-mm-radius anode 
exhibits nearly negligible variations. 

It is worth noting that the ratio (y/r,) can theoretically become a meaningful 
factor even for the 'normal-sized' anodes; extremely low-density plasmas could 
possibly generate values of y in the millimeter range. Thus, the notion of sheath 
characteristics being strictly a function of anode radius is valid only for the 


conditions being considered in this work. 
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Kplus, Cyl. Anode Comparison, Case | 
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Figure 45  K’ Profiles, a Three-Anode Comparison 
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(Kplusy, Cyl. Anode Comparison, Case I 
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Figure 46 (K*) Profiles, a Three-Anode Comparison 
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n-Profiles, Cyl. Anode Comparison, Case I 
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Figure 47 n,, Profiles, a Three-Anode Comparison 
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Kplus, Cyl. Anode Comparison, Case III 
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Figure 48K’ Profiles, a Three-Anode Comparison 
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(Kplusy, Cyl. Anode Comparison, Case II 
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Figure 49 (K’) Profiles, a Three-Anode Comparison 
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n-Profiles, Cyl. Anode Comparison, Case II 
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Figure 50 n,, Profiles, a Three-Anode Comparison 
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IV. CONCLUSIONS 


The manner in which one-dimensional planar and cylindncal anodes disturb 
collisional, isothermal plasmas has been explored in detail. Appropriate sheath and 
ambipolar regions have been shown to develop, for a given specific and reasonable 
electric field distribution, across a wide range of plasma types and conditions. 
Numerical solutions of the nondimensional governing equations have allowed specific 


and important observations to be made concerning: 
* the effect of temperature variation on plasmas in the anode region 


* traits that are characteristic of various plasma ‘families’, and traits that appear 
to be common to all collisional, low-temperature plasmas 


¢ the effect that differing anode radii have on the sheath for cylindrical anodes 


* the degree of difference between sheaths of planar and cylindrical anodes 
In addition, an analytical method has been presented that offers a simplified yet 
accurate algebraic solution for the lower-temperature plasmas. . 


The following briefly summarizes the specifics for each of these results. 


f 


Appropriate Sheath Existence: In each case considered, the {K’) curves showed 


that the net rate-of-charge-production decreases and vanishes toward the outer edge 
of the sheath. Moreover, each of the n-profiles clearly depicted sheath and ambipolar 
regions which were reasonably located with respect to the anode surface. 
Overall-Temperature_ Variation Effects: Results of this work indicate that at 
lower temperatures, charge production in the outer sheath is generic to the electric 
field distribution, and is thus independent of the details of ionization and 
recombination. Additionally, the n,, profiles, and thus the sheaths themselves, have 
been shown to be nearly unaffected by substantial changes in temperature. This 
somewhat unexpected result, which held for every plasma case considered, reduces 


any consequences caused by the isothermal-particle assumption (i.e., the assertion that 
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T. =T, =T,). It is interesting to note that the large magnitudes and physical meaning 
of n (the ratio of electrical energy to thermal energy in the sheath) actually predict 
such independence from temperature variations. 

Plasma ‘Family' Characteristics: All low-temperature, collisional plasmas can be 
partitioned into broad groupings based on their specific nondimensional ‘z’ or ‘q’ 
value (z denoting a planar anode problem, and q signifying cylindrical anodes). Each 
plasma ‘family’ thus exhibits its own characteristic electric field and sheath 
properties, even though its members may be widely diverse in composition and 
condition. 

Plasmas with large z (or q) values show significant charge production throughout 
the entire sheath, in contrast to the localized and decreased production rates depicted 
for the ‘small z/q' plasmas. This may indicate the fact that electric-field effects are 
more modest and thus localized in those plasmas with smaller values of z or q. In 
addition, the magnitude of the charged particles near the anode changes dramatically 
as ‘z/q' varies, especially in the sheath. The sheaths of ‘large z/q' plasmas exhibit 
charged-particle densities up to one order-of-magnitude smaller than those of the 
‘small z/q’ sheaths. More importantly, the sheath itself extends further from the 
anode surface as the value of 'z/q' decreases. In particular, the sheath of the ‘small 
z/q' plasmas stretches over ten times the distance occupied by the ‘large z/q' sheath. 

Cylindrical Anode Radius Effects: The sheath and the charge production rates 
have been shown to vary as the anode radius varies; however, the differences are 
almost negligible for all practical anodes and conditions. It appears that only 
extremely low-density plasmas, or cylindrical anodes whose radii are on the order of 


the sheath thickness itself, would exhibit significantly different characteristics. 
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Planar and Cylindrical Anode Differences: In an extension of the previous 
paragraph, the results of this work indicate that for the range of conditions explored, 
all practical cylindrical anodes can be treated to behave as planar anodes. 

Simple Analytical Approximations (Low Temperatures Only): Using the ‘Outer 
Solution’ technique, an algebraic solution can be used to generate highly accurate 
approximations to low-temperature anode problems. In addition, as long as 
derivatives of the presumed electric field distribution exist, the sheath profiles for all 


cases can be produced analytically to any desired accuracy. 


One suggestion for further work involves researching the effects that other 
electric field distributions would have on the anode problem. Such distributions may 
be derived from empirical data or from other ‘guesstimates', and may have three or 
more adjustable parameters. The most prominent reason for this suggestion 1s the 
need to corroborate the generic charge production rate observed at lower 
temperatures. | 

Another area requiring further study is the exploration of two-temperature 
plasmas and their anode sheaths. Results for such conditions could confirm the 


assertion of sheath insensitivity to temperature changes. 
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APPENDIX A 


Applicable FORTRAN Programs 
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Gatos PROGRAM IS BUILT AROUND A GRAFkKit 3.1 RUNGE-KUTTA 
meeoRIT THM PRODUCED BY SCO, INC. 


Pee PRECISION YA(0:10),YN(0:10),EK(0:4,0:10),Y¥(0:10),XA 
BevnbE PRECISION XP,XB,XM,H,HH,E,W, KPLUS,WPRM, KPLSPRM 

Pee PI,XL,XCHK,EINF,E0O,Z,ETA 

INTEGER NS,Pl 


ER INT * 

PeaNT *,’ FOURTH ORDER RUNGE-KUTTA SCHEME ’ 
PRaoNnT *,'’ FOR THE CARTESIAN COMPUTATION ' 
PaoiNrT *,’ OF W, KPLUS, & KPLSPRM’ 


Peete *,* INPUT VALUES FOR EINF, EO, Z, AND ETA’ 
READ *,EINF,E0O,Z,ETA 


IM=1 ! Number of equations 

meeee— 1.00 miniuwdumeondition for y 1 at x=XP. 
Yizew= 0 elie raimceondltron for y 2 at x=XP (if nec) 
PRINT * 

PRINT *,’ INTERVAL OF X FOR PRINTING ?’ 

READ *,PI 


Peer *,’ INPUT THE STEP SIZE (delta-x)’ 
READ *,H 


NS = NINT(PI/H) 


Beene *,’ MAXIMUM X TO STOP CALCULATION ?'° 


READ *,XL 
PRINT *,’ H= ',H 
Ea— © 

XP=0 

HH=H/2 

KPLUS=1.0 


print *,’ NS= , NS 
PeeiNT * 
LI = 0 ! Line no. initialization 
PRINT*, ’ LINE v7 a KPEUS” 
WRITE (*,98) LI,XP, KPLUS 
LI=LI+1 
DO N=1,NS 
XB=XP ! Old time 
XP=XP+H ! New time 
XM=XB+HH ! Midpoint time 
J=1 {! This part computes k°l. 
DO I=1,I1IM 
Val )j=YCl 
END DO 
XA=XB 
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CALL FUNCT( ERR ae, XA, BUN ECs eee 


J=2 ! This part computes k°2. 
DO i=l 7h 

YA(I)=Y(1I)+EK(1,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,E0O,Z,ETA) 
J=3 ! This part computes k7~3. 
DO I=1,I1M 

YA(I)=Y(1I)+EK(2,1)/2 
END DO 
XA=XKM 


CALL FUNCT(EK,J,YA,H,XA,EINF,E0O,Z,ETA) 


J=4 ! This part computes k4. 
DO I=1,I1M 
YA(I)=Y(1I)+EK(3,1) 
END DO 
XA=XP 
CALL FUNCT(EK,J,YA,H,XA,EINF,E0,Z,ETA) 
DO I=1,I1M ! 4-th order Runge-Kutta scheme 
Y(I)=Y(1I)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,1))/6 
END DO 
W=Y(1) 
You now have W( y/a ) -- to get K+( y/a ) you must 
multiply by { El y/a )°7 7E 00m, 
E=(EINF/EO) *DEXP((Z**2.)/((XP4+1.)%**2.)) 
KPLUS=W*E 
To get the derivative of W [ W’( y/a ) J, just plug the 
computed values of W back in the original lst-order ODE 
J=5 
YA(1)=W 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,E0,Z,ETA) 
WPRM=EK(5,1)/H 


NOW, to get d/dy [K+] from d/dy [ K+/(E/E0O) ] , must 
perform the following operation: 
KPLSPRM=(WPRM*(E**2.)+KPLUS*E*(-2.%*(Z**2.)/((XP+1.)**3.)))/E 


To keep from generating unplottable 50,000 point data files, 
the following will edit out data points depending on 
’where’ they occur 


XCHK=(XP/H) 

IF (XCHK .GT. 100000) GOTO 72 
IF (XCHK (GT. e000 C) ee to. 73 
IF (XCHK ;GT. 71000 GoTo 4 
IF (XCHK .GT.2 100) Geilo sys 
GOTO 44 


P1=P1+1 Ve 


IF (Pl .NE. 10000) GOTO 88 
P1=0 

GOTO 44 

P1=P1+1 

IF (Pl .NE. 1000) GOTO 88 
P1i=0 

GOTO 44 

Pl=P1+1 

imeeeeet .NE. 100) GOTO 88 
P1=0 

GOTO 44 

P1l=Pl1+1 

IF (Pl .NE. 10) GOTO 88 
P1=0 


WRITE (13,*) XP, KPLSPRM 
WRITE (12,*) XP, KPLUS 
WRITE (10,*) XP, W 


CONTINUE 


END DO 
WRITE (*,98) LI,XP, KPLUS 
Peper; lx, I2, F10.6, 2X, 1P4E16.8) 


eeeeeees .LT. XL) GOTO 28 


PRINT* 
PRINT*,’TYPE 1 TO CONTINUE, OR 0 TO STOP.’ 
READ *,K 

IF(K.EQ.1) GOTO 1 

PRINT* 

END 
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SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,EO,Z,ETA) ! DEFINES SET OF EQS 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),H,XA,PART1(0:4,0:10) 
DOUBLE PRECISION PART2(0:4,0:10) 
PART1(J,1)=DEXP((2.*2Z**2.)/((XA+1.)**2. 
PART2(J,1)=(DEXP((Z2**2.)/((XA4+1.)**2.)) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)-(ETA* ( 
EK(J,2)= the second ode, if nec. 
RETURN 

END 


))/((XA+1.)**3.) 
)* (ETA*EINF/EO) 
EINF/EOQ)**2.)*PART1(J,1))*H 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 


DOUBLE PRECISION YA(0:10),YN(0:10),EK(0:4,0: UOney (0 tOp ee 
DOUBLE PRECISION XP,XB,XM,H,HH,E,W,KPLUS,WPRM, KPLSPRM 
DOUBLE PRECISION NI,NE,NI1,NI2,N1I3,NI2A,NI3A,NE1 

REAL P1,XL,XCHK, EINE, BO, 2,2 le 

REAL EPSOQ,CHRG,K,TEMP,NINF,A 

INTEGER NS,P1l 


EPSO=8.853742E-12 
CHRG=1.602E-19 
K=1.38054E-23 


PRIN Te 

PRE) FOURTH ORDER RUNGE-KUTTA SCHEME ’ 
PRINT *,’ FOR THE CARTESIAN COMPUTATION '’ 

PRINTS OF NI AND NE’ 


PRINT *,’INPUT VALUES FOR EINF, EQ, Z, AND ETA’ 
READ *,EINF,E0,2Z,ETA 

PRINT *,’INPUT VALUES FOR TEMP. AND NINF’ 

READ *,TEMP,NINF 

PRINT *7 INPUL eva nUER Fr Ora. 


READ *,A 

IM=1 ! Number of equations 

Viel) -— 100 ! Initial condition for y 1 at x=XP. 

Vi en) ! Initial condition for y 2 at x=XP (if nec) 
PRINT * 

PRINT *,’INTERVAL OF X FOR PRINTING ?’ 

READ *,PI 


PRINT *,’INPUT THE STEP SIZE (delta-x)’ 
READ *,H 


NS = NINT(PI/H) 


PRINT *,’MAXIMUM X TO STOP CALCULATION ?' 
READ *,XL 


PRINT *,| )H-=Ve 


12 0) 

XP=0 
HH=H/2 
KPLUS=1.0 


LI = 0 ! Line no. initialization 
PRINT*, LINE y/a KPLUS’ 
WRITE (*,98) LI,XP, KPLUS 
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LI=LI+1 


DO N=1,NS 
XB=XP ! Old time 
XP=XP+H ! New time 
XM=XB+HH ! Midpoint time 
J=1 ! This part computes kl. 
DO I=1,IM 
YA(I)=Y(I) 
END DO 
XA=XB 
CALL FUNCT(EK,J,YA,H,XA,EINF,E0,2Z,ETA) 
J=2 ! This part computes k°2. 
DO I=1,I1IM 
YA(I)=Y(1I)+EK(1,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,E0O,Z,ETA) 
J=3 ! This part computes k°3. 
DO I=1,IM 
YA(I)=Y(1I)+EK(2,1)/2 
END DO 
XA=XM 


Sabi FUNCT(EK,J,YA,H,XA,EINF,E0,Z,ETA) 


J=4 ! This part computes k°4. 
DO I=1,IM 
YA(I)=Y(1I)+EK(3,1I) 
END DO 
XA=XP 
CALL FUNCT(EK,J,YA,H,XA,EINF,E0,Z,ETA) 
DO I=1,IM ! 4-th order Runge-Kutta scheme 
Y(I)=Y¥(1I)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,1))/6 
END DO 
W=Y(1) 
You now have W( y/a ) -- to get K+( y/a ) you must 
multiply by te E( y7a } 7 E(O) } 
E=(EINF/EO) *DEXP( (Z**2.)/((XP+1.)**2.)) 
KPLUS=W*E 
To get the derivative of W [ W’( y/a ) J, just plug the 
computed values of W back in the original lst-order ODE 
J=5 
YA(1)=W 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,E0,Z,ETA) 
WPRM=EK(5,1)/H 


NOW, to get d/dy [K+] from d/dy [{ K+/(E/E0O) ] , must 
perform the following operation: 
KPLSPRM=(WPRM*(E**2.)+KPLUS*E*(-2.*(Z**2.)/((XP+1.)**3.)))/E 
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AND NOW, to generate the non-dimensionalized n-curves 
ne( y/fa ) and nz( ya ) 
just compute these equations: 


NI1=(KPLUS/E) *(EINF/EO) 
NI2A=(EPSO*K*TEMP) /( (CHRG**2.) *NINF*(A**2.)) 
NI2=NI2A*((6.*(Z2**2.)/((XP4+1.)**4.))4(4.%(Z2**4.)/((XP41.)**6.))) 
NI3A=((EPSO*EQ) /(CHRG*NINF*A) ) *E 
NI3=NI3A*(-2.*(Z**2.)/((XP41.)**3.)) 

NI=.5*(NI14+NI2+NI3) 
NE1=((EQO*EPSO*(Z**2.))/(CHRG*NINF*A) ) *E*(2./((XP4+1.)**3.)) 

NE=NI+NE1 


To prevent every data point from being written to the file 
(resulting in unplottable 50,000 pt. data files), the following 
edits out a percentage of the data depending on ‘’where’ it 
was generated. 


XCHK=(XP/H) 

IF (XCHK .GT. 100000) GOTO 72 
IF (XCHK .GT. 10000) GOTO 73 
IF (XCHK .GT. 1000) GOTO 74 
IF (XCHK .GT. 100) GOTO 75 
GOTO 44 

P1=P1+1 

IF (Pl .NE. 10000) GOTO 88 
P1=0 

GOTO 44 

P1=P1+1 

IF (Pl .NE. 1000) GOTO 88 
Pl=0 

GOTO 44 

P1=P1+1 

IF (Pl .NE. 100) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 10) GOTO 88 
Pl=0 


WRITE (007 Xe ee 
WRITE (11,*) XP, NE 
CONTINUE 
END DO 
WRITE (*,98) LI,XP, KPLUS 


FORMAT(IX, i2, FIO .6,5 2k, 1 PdElowu® 


IF (XP .LT. XL) Gores 


PRINT 

PRINT*,’TYPE 1 TO CONTINUE, OR 0 TO STOP.” 
READ *,Kl 

TP(K1. EO] CO Lom 

Jer agit Nib os 
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END 
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SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,E0O,Z,ETA) ! DEFINES SET OF EQS 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),H,XA,PART1(0:4,0:10) 

DOUBLE PRECISION PART2(0:4,0:10) 
PART1(J,1)=DEXP((2.*Z**2.)/((XA+1.)**2.))/((XA+1.)**3.) 
PART2(J,1)=(DEXP((Z**2.)/((XA+1.)**2.)))*(ETA*EINF/EO) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)-—(ETA*(EINF/EO)**2.)*PART1(J,1))*H 
EK(J,2)= the second ode, if nec. 

RETURN 

END 


We, 


THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED=EY SCO, INGA 


DOUBLE PRECISION YA(0:10),YN(0:1OQ)  EK(0:4,0: 1000; 10h 
DOUBLE PRECISION XP,XB,XM,H,HH,E,W, KPLUS,WPRM, KPLSPRM 

REAL PI,XL,XCHK,EINF,ERO,Q,ETA 

INTEGER NS,Pl 


PRINT * 

PRINTS 7 FOURTH ORDER RUNGE-KUTTA SCHEME ’ 
BRUNT i FOR THE CYLINDRICAL COMPUTATION ~ 
PRON OF W, KPLUS, & KPLSPRM’ 

PRINT =, ( NORMAL-SIZED ANODE: r0 = 1 cm )’ 


PRINT *,’INPUT VALUES FOR EINF, ErQ, Q, AND ETA’ 
READ *,EINF,ERO,Q,ETA 


IM=1 ! Number of equations 

Y(i) = 1200 ! Initial condition for wl at r° =XP. 

YaC7 ea) ! Initial condition for w2Zvat, =x ee nee) 
PRINT * 

PRINT *,’INTERVAL OF r” FOR PRINTING ?’ 

READ *,PI 


PRINT *,’INPUT THE STEP SIZE (delta-r )’ 
READ *,H 


NS = NINT(PI/H) 


PRINT *,’MAXIMUM r”~ TO STOP CALCULATION ?’ 
READ *,XL 


PRINT *,’ dr~ = ',H 
Pi=0 
XP=0 


HH=H/2. 
KPLUS=1.0 


print * ,; No=a Ns 


PRINT * 
Jl 5), ! Line no. initialization 
PRIN LINE ri KPLUS’ 


WRITE (*7958) LrSXe eke ru. 


LI=LI+1 
DO N=1,NS 
XB=XP ! Old time 
XP=XP+H ! New time 
XM=XB+HH 1 Midpoint time 
Sa Runge-Kutta Scheme---~-~~-~------------ 
J=1 Thins swan eee 
DO I=1,IM 
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weet I )=Y(1) 
END DO 
XA=XB 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 


J=2 ww ris 1S part. 2" 
DO I=1,1M 

YA(1I)=Y(1I)+EK(1,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 
J=3 (ernas iS part 3. 
DO I=1,IM 

me )=Y(1)+EK(2,1)/2 
END DO 
XA=XM 

CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 
J=4 ! This 1s part 4. 
DO I=1,I1M 

YA(I)=Y(1I)+EK(3,1) 
END DO 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 


mo, t=1,IM ! 4-th order Runge-Kutta scheme 
Y(I)=Y¥(1I)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,1))/6 
END DO 
W=Y(1) 
You now have W( r” ) -- to get K+( xr” ) you must 
multiply by f,E( ©” ) / Ee 4 
E=(EINF/ERO ) *DEXP((Q**2.)/((XP+1.)**2.)) 
KPLUS=W*E 
To get the derivative of W [ W’( xr” ) J, just plug the 
computed values of W back in the original lst-order ODE 
J=5 
ee) =W 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 
WPRM=EK(5,1)/H 


NOW, to get d/dr”™ [K+] from d/dr”™ [ K+/(E/Er0) J] , must 
perform the following operation: 
KPLSPRM=(WPRM*(BE**2.)4+KPLUS*E*(-2.%*(Q**2.)/((XP+1.)**3.)))/E 


To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on ‘where’ they occur 


meecHK=(XP/H) 
meet XCHK .GT. 100000) GOTO 72 


mee (XCHK .GT. 10000) GOTO 73 me 


IF (XCHK .GT. TOW wGeToO 72 
IF (XCHK .GT.” TOUpeiGorO 75 
GOTO 44 


FZ P1l=P1+1 
IF (Pl .NE. 10000) GOTO 88 
P1=0 
GOTO 44 

13 Pl=P1+1 
IF (Pl .NE. 1000) GOTO 88 
P1=0 
GOTO 44 

We P1l=P1+1 
IF (Pl .NE. 100) GOTO 88 
P1=0 
GOTO 44 

Tes P1l=P1+1 
IF (Pl .NE. 10) GOTO 88 
Pi 


44 WRITE (137%) XEe KPLSPRM 
WRITE (1254 )exe, KPLUS 
WRITE (14,*) XP, W 


88 CONTINUE 


END DO 


WRITE (*,98) LI,XP, KPLUS 
98 FORMAT( 1X, 12, F110 56,5 ZX I P4E lowe 


IF (XP .LT. XL) Gore Ze 


200 PRINT* 
PRINT*, ‘TYPE 1 TO CONEINUE, OR OG) TOmsiere 
READ *,Kl 
PreK. EO.1) GoTor! 
PRINT* 
END 
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SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) ! DEFINES SET OF EQS 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),H,XA,PART1(0:4,0:10) 
DOUBLE PRECISION PART2(0:4,0:10) 
PART1(J,1)=ETA*(DEXP((2.*Q**2.)/((XA+1.)**2.))/((XA+1.)**3.)) 
PART2(J,1)=(DEXP((Q**2.)/((XA+1.)**2.)))*(ETA*EINF/ERO) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)-((EINF/ERO)**2.)*PART1(J,1))*H 

C EK(J,2)= the second ode, if nec. 
RETURN 
END 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 


Pemeoe PRECISION YA(0:10),YN(0:10),EK(0:4,0:10),Y(0:10),XA 
DOUBLE PRECISION XP,XB,XAM,H,HH,E,W,KPLUS,WPRM, KPLSPRM 
DOUBLE PRECISION NI,NE,NI1,N1I2,N13,NE1 

meee! ,CC2,€C3,82,83,S82A,S2B,S82C 

REAL PI,XL,XCHK,EINF,ERO,Q,ETA 

REAL EPSO,CHRG,K,TEMP,NINF,B,RO 

INTEGER NS,Pl 


EPSO=8 .853742E-12 
CHRG=1 .602E-19 
K=1.38054E-23 


RO=.01 

PFRIENT * 

Ereeua *,' FOURTH ORDER RUNGE-KUTTA SCHEME ' 
Perna *, ' FOR THE CYLINDRICAL COMPUTATION ’ 
ERINT *,' OF NE AND NI’ 

FeuNT *,’ ( NORMAL-SIZED ANODE: EU =— sem )? 


Prete,’ INPUT VALUES FOR EINF, Er0O, Q, AND ETA’ 
READ *,EINF,EROQ,Q,ETA 

PRINT *,’INPUT VALUES FOR TEMP. AND NINF’ 

READ *,TEMP,NINF 

PRINT *,’INPUT VALUE FOR b’ 


READ *,B 

IM=1 ! Number of equations 

Mel) = 1.00 minveral Condition for wl at rc” =XP. 

vaee,) = 0 [emer cial condition for wZ2 at r =XP (if nec) 
PRINT * 

PRINT *,’INTERVAL OF r’ FOR PRINTING ?° 

READ *,PI 


PRT *,’ INPUT THE STEP SIZE (delta-r’)’ 
READ *,H 


NS = NINT(PI/H) 


PRINT *,’MAXIMUM r”~ TO STOP CALCULATION ?' 
READ *,XL 


Poeun *,’ dr’ = ',H 


P1=0 
XP=0 
HH=H/2. 
KPLUS=1. 


Cie 0 ! Line no. initialization 


FRINT*, ’LINE E KEES 


WRITE (*,98) LI,XP, KPLUS 81 
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Li= List 


DO N=1,NS 
XB=XP ! Old time 
XP=XP+H ! New time 
XM=XB+HH ! Midpoint time 
—--------- Runge-Kutta Scheme--—-—---—------------- 
J=1 ! This is part l. 
DO I=1,IM 
YA(I)=Y(1I) 
END DO 
XA=XB 
CALL FUNCT(EK,J,YA,H,XA, EINF,EROQ,Q,ETA) 
J=2 ! This is part 2. 
DO I=1,IM 
YA(I)=Y(1I)+EK(1,1)/2 
END DO 
XA=KM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 
J=3 i THis .1S.pamt 3. 
DO I=1,IM 
YA(I)=Y(I)+EK(2,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q, ETA) 
J=4 ! This is part 4. 
DO I=1,IM 
YA(I)=Y(1I)+EK(3,1) 
END DO 
XA=XP 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) 
DOm l= een ! 4-th order Runge-Kutta scheme 
Y(I)=Y(1I)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,I1))/6 
END DO 
W=Y(1) 
You now have W( r” ) -- to get K+( r™~ ) you must 
multiply by { Et “2 JE CO) 
E=(EINF/ERO) *DEXP((Q**2.)/((XP+1.)%**2.)) 
KPLUS=W*E 
To get the derivative of W [ W’( r~ ) ], just plug the 
computed values of W back in the original l1st-order ODE 
J=5 
YA(1)=W 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA, BINF, ERO, Q, ETA) 


Newe to get d/dr [K+] from d/dr° [ K+/(E/Er0) ] , must 
perform the following operation: 
KPLSPRM=(WPRM*(E**2.)+KPLUS*E*(-2.*(Q**2.)/((XP+1.)**3.)))/E 


AND NOW, to generate the nondimensionalized n-curves 
ne( rc’ ) and ni( r” ) 
compute the following equations 


NI1=(KPLUS/E) *(EINF/ERO) 
CC2=(K*TEMP*EPSO) /( (CHRG**2.)*NINF*(B**2.)) 
(6. *(0**%2.))/((XP4+1.)**4. ) 
S2B=(4.*(Q**4.))/((XP41.)**6.) 
S2C=(1/(XP+(RO/B)))*((-2.*(Q**2.))/((XP+1.)**3.)) 
S2=S2A+S2B+S2C 


NI2=CC2*S2 
CC3=(EPSO*ERO)/( CHRG*NINF*B) 
S3=((-2.*(Q**2.))/((XP41.)**3.))4+(1./(XP+(RO/B) ) ) 


NI3=CC3*E*S 3 
NI=.5*(NI1+NI2+NI3) 


CC1l=(EPSO0*ERO) /( CHRG*NINF*B) 
NE1L=CC1*E*((-2.*(Q**2.))/((XP4+1.)**3.)4+(1./(XP+(RO/B) ) ) ) 


NE=NI-NE1 


To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on ‘/where’ they occur 


XCHK=(XP/H) 

IF (XCHK .GT. 100000) GoTo 72 
IF (XCHK .GT. 10000) GoTO 73 
IF (XCHK .GT. 1000) GOTO 74 
me (XCHK .GT. 100) GOTO 75 
GOTO 44 


P1l=P1+1 

IF (Pl .NE. 10000) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 1000) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

feet Pl .NE. 100) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (P1 .NE. 10) GOTO 88 
P1=0 


WRITE (10,*) XP, NI 
WRITE (11,*) XP, NE 


CONTINUE 
83 


END DO 


WRITE (*,98) LI,XP, KPLUS 
98 FORMAT(1X, I2, FIEUS6R2x, Sedcloaee 


IF (XP -LT. XL)eGere 23 


200 PRINT* 
PRINT*, ‘TYPE 1 TO CONTINUE, JORMO.10n foc. 
READ *,K1l 
TECK EO. eGOTer! 
PRIN 
END 
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SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA) ! DEFINES SET OF EQS 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),H,XA,PART1(0:4,0:10) 
DOUBLE PRECISION PART2(0:4,0:10) 
PART1(J,1)=ETA*(DEXP((2.*0**2.)/((XA+1.)**2.))/((XA+1.)**3.)) 
PART2(J,1)=(DEXP((Q**2.)/((XA+1.)**2.)))*(ETA*EINF/ERO) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)-((EINF/ERO)**2.)*PART1(J,1))*H 

Gc EK(J,2)= the second ode, if nec. 
RETURN 
END 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 


Deuvene PRECISION YA(0:10),YN(0:10),EK(0:4,0:10),Y¥(0:10),XA 
DOUBLE PRECISION XP,XB,XM,H,HH,E,W,KPLUS,WPRM, KPLSPRM 

REAL PI,XL,XCHK,EINF,ERO,Q,ETA,RO,G,B,EPS0,CHRG 

INTEGER NS,P1l 


PRINT * 

Pein? *, ’ FOURTH ORDER RUNGE-KUTTA SCHEME ' 
PeiINT *,' FOR THE CYLINDRICAL COMPUTATION '’ 
Pot *, ’ OF W, KPLUS, & KPLSPRM’ 

Peat *,’ ( WIRE-THIN ANODE: eae in 


PRINT *,’INPUT VALUES FOR EINF, Er0Q, Q, AND ETA’ 
READ *,EINF,ERO,Q,ETA 

PRINT *,’INPUT VALUES FOR GAMMA AND B’ 

READ *,G,B 


R0=.0001 
EPSO=8 .853742E-12 
CHRG=1 .602E-19 


PM=1 ! Number of equations 

Yil) = 1.00 minitwval condition ror wl ae r =XP. 

viene = 0 elnmitialweeoendition for w2 at fr =XP (if nec) 
PRINT * 

PRINT *,’INTERVAL OF r” FOR PRINTING ?’ 

READ *,PI 


Poe *,’ INPUT THE STEP SIZE (delta-r” )’ 
READ *,H 


NS = NINT(PI/H) 


PRINT *,’MAXIMUM r~ TO STOP CALCULATION ?’ 
READ *,XL 


PRINT *,’ dr~ = ’,H 


2 10 
XP=0 
HH=H/2. 
KPLUS=1. 


print *,’ NS= ’, NS 
PRINT * 


ro = QO ' Line no. initialization 
PRINT*, ’LINE Tees KPLUS’” 
WRITE (*,98) LI,XP, KPLUS 


LI=LI+1 

DO N=1,NS 
XB=XP Old stime 
XP=XP+H ! New time 
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XM=XB+HH ! Midpoint time 


J=1 iO this 21s gpartes. 
DO I=1,I1M 
YA(I)=Y(I1) 
END DO 
XA=XB 
CALL FUNCT( EK QJ YA ,H,23A, EDNE ERO OETA Gee. 
J=2 PSTn1iS#isS pate z- 
DO I=1,I1IM 
YA(I)=Y(I)+EK(1,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
J=3 i This as Part 3. 
DO I=1,IM 
YA(I)=Y(1)+EK(2,1)/2 
END DO 
XA=XKM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
J=4 ! This 2s pares 
DO I=1,I1M 
YA(I)=Y(1)+EK(3,1) 
END DO 
XA=XP 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
DO yl—ia ! 4-th order Runge-Kutta scheme 
Y(I)=Y(1)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,1))/6 
END DO 
W=Y(1) 
You now have W( r~ ) -- to get K+( r~ ) you must 
multiply bx, { E(t sc” ) 7 Ee 


F=(EINF/ERO) *DEXP((Q**2.)/((XP+1.)**2.)) 
KPLUS=W*E 


To get the derivative of W [ W’( r” ) J, just plug the 
computed values of W back in the original 1st-order ODE 


=> 

YA(1)=W 

XA=XP 

CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
WPRM=EK(5,1)/H 


NOW, to get d/dr™~ [K+] from d/dr~ [ K+/(E/Er0) ] , must 
perform the following operation: 
KPLSPRM=(WPRM*(E**2.)+KPLUS*E*(-2.*(Q**2.)/((XP+1.)**3.)))/E 


To keep from generating unplottable 50,000 point data files, 
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the following edits out some of the data points depending 
on ‘where’ they occur 


XCHK=(XP/H) 

Maen CHK .GT. 100000) GOTO 72 
MmemeexCHK .GT. 10000) GOTO 73 
meecxCHK .GT. 1000) GOTO 74 
Mee xCHK .GT. 100) GOTO 75 
GOTO 44 


P1=P1+1 

IF (Pl .NE. 10000) GOTO 88 
Ea =( 

GOTO 44 

P1l=P1+1 

meetel .NE. 1000) GOTO 88 
Pi1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 100) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 10) GOTO 88 
P1=0 


WRITE (13,*) XP, KPLSPRM 
WRITE (12,*) XP, KPLUS 
WRITE (14,*) XP, W 


CONTINUE 


END DO 


WRITE (*,98) LI,XP, KPLUS 
Boemmar( ix, I2, F10.6, 2X, 1P4E16.8) 


Maeeene .LT. XL) GOTO 28 


PRINT* 

eee TYPE 1 TO CONTINUE, OR O TO STOP.’ 
READ *,K1l 

tah. EQ.1) GOTO 1 

Een T * 

END 


EK KKKKKKEKKEKEKKEKEKEKEKKKKKKKKKK KKK KKKKKKKKE 


SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) ! DEFINES ODEs 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),PART2(0:4,0:10) 
DOUBLE PRECISION PART3(0:4,0:10),PART4(0:4,0:10),H,XA,Cl1l 
REAL RO 

ro=- 0001 

C2=ETA*G*((EINF/ERO) **2. ) 
PART4(J,1)=C2*(DEXP((2.*(Q**2.))/((XA+1.)**2.)))/(XA*B+R0) 
C1=ETA*((G/RO)+1.)*( (EINF/ERO)**2. ) 
PART3(J,1)=C1*(DEXP((2.*Q**2.)/((XA41.)**2.))/((XAt1.)**3.)) 
PART2(J,1)=(DEXP((Q**2.)/((XA+1.)**2.)))*(ETA*EINF/ERO) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)-PART3(J,1)+PART4(J,1))*H 


87 


é 


EK(J,2)= 
RETURN 
END 


the second ode, 


if nec. 


88 


THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 


Pewee PRECISION YA(0:10),YN(0:10),EK(0"4,0:10),Y(0:10),XA 
DOUBLE PRECISION XP,XB,XM,H,HH,E,W,KPLUS,WPRM,KPLSPRM 
DOUBLE PRECISION NI,NE,NI1,NI2,N13,NEl1 

pee 1,CC2,CC3,82,83,82A,S2B8,82C,S82D 

Ree ri ,XL,XCHK,EINF,ERO,Q,ETA 

REAL EPSO,CHRG,K,TEMP,NINF,B,G,RO 

INTEGER NS,P1l 


EPSO=8 .853742E-12 
CHRG=1.602E-19 
K=1.38054E-23 


RO=.0001 

PRINT * 

PRs *, ’ FOURTH ORDER RUNGE-KUTTA SCHEME ’ 

Pane *,' FOR THE CYLINDRICAL COMPUTATION '’ 

Pein *,' OF NE AND NI’ 

Eee NT *, ’ ( WIRE-THIN ANODE: r0 = 0.1 mm )’ 


PRINT *,’INPUT VALUES FOR EINF, Er0O, Q, AND ETA’ 
Peeb *,EINF,ERO,Q,ETA 

Peon *,’INPUT VALUES FOR TEMP. AND NINF’ 

READ *,TEMP,NINF 

PRINT *,’INPUT VALUES FOR GAMMA AND b’ 

READ *,G,B 


IM=1 ! Number of equations 

weet) = 1.00 minitial cenmdition Lor wl at rr =XP. 

ec) = O Meinitral Condition For wZ at r =XP (if nec) 
PRINT * 

PRINT *,’INTERVAL OF r”~ FOR PRINTING ?’ 

READ *,PI 


PRINT *,’INPUT THE STEP SIZE (delta-r’ )’ 
READ *,H 


NS = NINT(PI/H) 


PRINT *,’MAXIMUM r” TO STOP CALCULATION ?’ 
READ *,XL 


Peans *,'’ dr” = ’',4H 


P1=0 
XP=0 
HH=H/2. 
KPLUS=1. 


print *,’ NS= ’, NS 
PRINT * 


ho 0 PeLine uno ....1nicaialization 
PRINT*, ’LINE a KPLUS’ 
WRITE (*,98) LI,XP, KPLUS 
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28 


LI=LI+1 


DO N=1,NS 
XB=XP ! Old time 
XP=XP+H ! New time 
XM=XB+HH ! Midpoint time 


---------- Runge-Kutta Scheme------------------- 


J=1 ! This is part l. 
DO I=1,IM 
YA(I)=Y(I) 
END DO 
XA=XB 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
J=2 PeCthRiS 1S epee 2. 
DO I=1,I1M 
YA(I)=Y(1)+EK(1,1)/2 
END DO 
XA=XM 
CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
J=3 "SHhLS 1S hear tase 
DO I=1,I1IM 
YA(1 )=Y( TP)4EK G22 2 
END DO 
XA=XM 


CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 


J=4 i This a&s part <4. 
DO I=1,I1IM 
YA(I)=Y(1I)+EK(3,1) 
END DO 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 


DO I=1,IM ! 4-th order Runge-Kutta scheme 
Y(I)=Y¥(1)+(EK(1,1)+EK(2,1)*2+EK(3,1)*2+EK(4,1))/6 
END DO 
W=Y(1) 
You now have W( rr” ) -- to get K+( rr” ) you must 
multiply by { E( 4c  )} 7 EGO 
E=(EINF/ERO) *DEXP( (Q**2.)/((XP+1.)**2.)) 
KPLUS=W*E 
To get the derivative of W [W’( xr” ) J, just plug the 
computed values of W back in the original ist-order ODE 
J=5 
YA(1)=W 
XA=XP 


CALL FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 
WPRM=EK(5,1)/H 
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Hewemco get d/dr [K+] from d/dr™ [{ K+/(E/Er0) J] , must 
perform the following operation: 
KPLSPRM=(WPRM*(E**2.)+KPLUS*E*(-2.*(Q**2.)/((XP+1.)**3.)))/E 


AND NOW, to generate the nondimensionalized n-curves 
ne( r” ) and ni( or” ) 
compute the following equations 


NI1=(KPLUS/E) *(EINF/ERO) 
CC2=(K*TEMP*EPSO) /( (CHRG**2.)*NINF*(B**2. ) ) 
S2A=(6.*(O**2.))/((XP41.)**4.) 
S2B=(4.*(Q**4.))/((XP+1.)**6.) 
S2C=(1./(XP+(RO/B) ))*((-2.*(Q**2.))/((XP+1.)**3.)) 
S2D=(1./((XP+(RO/B) )**2.)) 
S2=S2A+S2B+S2C-S2D 
NI2=CC2*S2 
CC3=(EPS0O*ERO) /( CHRG*NINF*B) 
§3=((-2.*(Q**2.))/((XP41.)**3.))+(1./(XP+(RO/B) ) ) 
NI3=CC3*E*S3 


NI=.5*(NI1+NI2+NI3) 


CC1l=(EPS0O*ERO) /( CHRG*NINF*B) 
NEI=CC1*E*((-2.*(Q**2.))/((XP41.)**3.)4(1./(XP+(RO/B) ))) 


NE=NI-NE1 


To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on ‘where’ they occur 


XCHK=(XP/H) 

Mme xcCHK .GT. 100000) GOTO 72 
mex cHK  .GT. 10000) GOTO 73 
mee (XCHK .GT. 1000) GOTO 74 
MeeeencHK .GT. 100) GOTO 75 
GOTO 44 


P1=P1+1 

IF (Pl .NE. 10000) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 1000) GOTO 88 
P1=0 

GOTO 44 

P1=P14+1 

IF (Pl .NE. 100) GOTO 88 
P1=0 

GOTO 44 

P1l=P1+1 

IF (Pl .NE. 10) GOTO 88 
Pl=0 


WRITE (10,*) XP, NI 
WRITE (11,*) XP, NE 


9 : 


88 CONTINUE 


END DO 


WRITE (*,98) L1I,XP, KELUS 
98 FORMAT(1X, 12, FROV67 2k) bese oer 


IF (XP .LT. XL} GOTO =ZzsB 


200 Praia 
PRINT*,’TYPE 1 TO CONTINUE, OR SUeoOeS TOL 
READ *,K1l 
Il Fhe BO. 1) GOTOa! 
PRINT* 
END 


CHEEK KKKKKEKKEKKKKEKKEKAKKAEKKEKAEKKKEKKAEKKEKKKKKEKKAEKKK 


SUBROUTINE FUNCT(EK,J,YA,H,XA,EINF,ERO,Q,ETA,G,B) ! DEFINES ODES 
DOUBLE PRECISION EK(0:4,0:10),YA(0:10),H,XA,PART2(0:4,0:10) 
DOUBLE PRECISION PART3(0:4,0:10),PART4(0:4,0:10) 
REAL RO,KK4,KK3 
RO=.0001 
KK4=ETA*G*( (EINF/ERO) **2. ) 
PART4(J,1)=KK4*(DEXP((2.*(Q**2.))/((XA+1.)**2.)))/(XA*B+RO0) 
KK3=ETA*((G/RO)+1.)*( (EINF/ERO) **2. ) 
PART3(J,1)=KK3*(DEXP((2.*Q**2.)/((XA¢1.)**2.))/( (XA+1.)**3.)) 
PART2(J,1)=(DEXP((Q**2.)/((XA+1.)**2.)))*( ETA*EINF/ERO) 
EK(J,1)=(2.*ETA-PART2(J,1)*YA(1)—-PART3(J,1)+PART4(J,1))*H 

C EK(J,2)= the second ode, if nec. 
RETURN 
END 
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